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ABSTRACT. 

A  numerical  model  which  utilizes  the  isobaric  vorticity  equation 
is  developed  and  applied  to  cumulus-scale  data.   The  model,  together 
with  a  modified  version  of  the  cumulus  convection  model  of  Weinstein 
and  Davis,  is  applied  to  data  obtained  from  the  National  Severe  Storms 
Laboratory  in  Norman,  Oklahoma.   The  calculations  yield  real  time  pre- 
dictions for  height,  temperature  and  relative  humidity  at  seven  pres- 
sure levels,  which  are  then  used  as  input  to  the  cumulus  convection 
model  to  obtain  vertical  profiles  of  various  parameters  at  specified 

grid  points . 

Some  results  of  the  calculations  are  presented  along  with  sug- 
gestions for  further  testing  and  improvement.   The  results  indicate 
that  further  modifications  to  the  approach  used  are  necessary  in  order 
to  provide  more  accurate  forecasts.   Values  of  the  individual  terms 
in  the  vorticity  equation  are  presented  as  computed  from  the  observed 
mesoscale  data. 


LIBRARY 

NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY.    CALIF.    93940 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  939434101 


TABLE   OF  CONTENTS 


Page 

Abstract  2 

List    of  Tables  5 

List    of   Illustrations  7 

Table    of   Symbols    and  Abbreviations  9 

Acknowledgement  11 

1.  Introduction  13 

2.  Modifications  to  the  W.D.  Model  17 

3.  Mathematical  Development  of  a  Forecast  Model  23 

A.  The  vorticity  equation  23 

B.  Vertical  motion  computations  27 

C.  Temperature  and  relative  humidity  forecast  scheme     29 

D.  Boundary  conditions  30 

E.  Time  steps  32 

4.  Initialization  of  the  Model  34 

5.  Results  of  the  Computations  42 

6.  Suggestions  for  Future  Developments  and  Improvements  53 

7 .  Summary  57 
References  58 
Appendices : 

A.  The  NSSL  Radiosonde  Network  59 

B.  On  the  Ellipticity  Requirement  of  the  Balance         60 
Equation 


C.   Listing  of  Computer  Program 


62 


J00H083TM> 


LIST   OF  TABLES 

Table  Page 

1  Table    of   Symbols    and   Abbreviations  9 

2  Minimum-Maximum  values    of    individual    terms  43 
in  height    prediction   equation    for    1700Z, 

30  May    1967 

3  Minimum-Maximum  values  of  individual  terms         44 
in  height  prediction  equation  for  2130Z, 

30  May  1967 

4  Sample  of  values  of  the  individual  terms  45 
in  height  prediction  equation 

5  Maximum  upward  and  downward  values  of  47 
vertical  motion  at  1700Z  and  2130Z 

6  Average  absolute  errors  of  height  fields  49 
(meters)  for  persistence  and  model  at 

1700Z  and  2130Z  for  all  seven  isobaric 
leve Is 


LIST   OF   ILLUSTRATIONS 

Figure  Page 

1  Non-geos trophic    nature    of   data  14 


2  Variation  of    (C      -C    )   with   temperature  20 

pv  w  r 

3  "Regions  of  Influence"  for  initialization  35 

4a  Hand  analysis  for  950  mb  at  1700Z  on  30  May        38 

1967 

4b  Computer  objective  analysis  for  950  mb  at  39 

1700Z  on  30  May  1967 

5a  Hand  analysis  for  500  mb  at  1700Z  on  30  May        40 

1967 

5b  Computer  objective  analysis  for  500  mb  at  41 

1700Z  on  30  May  1967 

6  Sample  of  output  from  W.D.  model  51 


TABLE  OF  SYMBOLS  AND  ABBREVIATIONS 


Threshold  value  below  which  conversion  from  cloud  water  to 

-3 
hydrometeor  water  does  not  take  place  (gm  m   ) . 


Threshold  value  below  which  conversion  from  cloud  water  to 
hydrometeor  water  does  not  take  place  (gm  Kg   ) . 


C  Specific  heat  of  water  vapor  at  constant  pressure 


pv 


vc 


.      -1  o-l. 
(cal  gm    K   )  . 


C  Specific  heat  of  liquid  water  (cal  gm    K   ). 

w  r 

e  Saturation  vapor  pressure  in  millibars. 
S  tr         r 

f  Coriolis  parameter. 

f  Mean  Coriolis  parameter  for  region  studied  =  f0c' 
m  jj 

g  Acceleration  of  gravity  (m  sec   ). 

J  Jacobian. 

Lf  Latent   heat    of    fusion    (cal  gm      ). 

L  Latent   heat    of   sublimation    (cal  gm      ). 
s 

L  Latent    heat    of   vaporization    (cal  gm      ). 

-3 

m  Kessler's    saturation  vapor   density    (gm  m      ). 

_3 

M  Kessler's    precipitation   content    (gm  m      ). 

p  Pressure    in  millibars    (mb) . 

Aq  Saturation  mixing    ratio   over  water  minus    saturation  mixing 
ratio   over    ice    (gm  gm      ). 

Q  Cloud  water   content    (gm  gm      ). 


Hydrometeor  water   content    (gm  gm      ) . 

Gas    constant    for   dry  air  ~  0.28704    (joules    gm  K      ) 

Relative  humidity  (fo)  . 


TABLE  OF  SYMBOLS  AND  ABBREVIATIONS,  CONTINUED 


t  Time . 

T  Temperature. 

u  Zonal  component  of  horizontal  velocity  (m  sec   ). 

v  Meridional  component  of  horizontal  velocity  (m  sec   ). 

V  Terminal  velocity  of  hydrometeors  (m  sec   ). 
\V  Horizontal  wind  vector. 

\V  Divergent  component  of  the  horizontal  wind  vector. 

A. 

\V.         Rotational  component  of  the  horizontal  wind  vector. 

w  Vertical  component  of  velocity  in  (x,  y,  z,  t)  coordinate 

-In 

system  (m  sec   ) . 

z  Height  of  isobaric  surface  in  meters. 

C  Ratio  of  gas  constants  ~  .622. 

Q  Vertical  component  of  relative  vorticity  (sec   ). 

_3 
p  Air  density  (gm  m   ). 

$  Geopotential  =  gz . 

2  -1 

i|f  Stream   function    (m     sec      ). 

to  Vertical   component    of  velocity   in    (x,    y,    p,    t)    coordinate 

system  =  -j~   (mb   hr        or  mb   sec      ). 

S  ci 

V  "Del"   operator   =  li-r 1-  rKr-. 


V 


2  .      .       .  d2  d2 

Laplacian   operator   -  — ~  H r- 

dx    dy 


2 

V  Finite  difference  form  of  the  Laplacian  operator. 

X  Average  value  of  any  parameter  X  between  two  pressure  levels 

(X)  Average  value  of  any  parameter  X  with  respect  to  temperature 


10 


ACKNOWLEDGEMENT 

The  author  wishes  to  thank  Dr.  R.  L.  Alberty  for  his  guidance  and 
encouragement  throughout  the  research  and  during  the  writing  of  this 
thesis,  Dr.  G.  J.  Haltiner  and  Dr.  J.  D.  Mahlman  for  their  invaluable 
suggestions  and  advice,  and  my  wife  and  children  for  their  patience 
throughout  the  project. 


11 


1.       INTRODUCTION 

The   National   Severe    Storms    Laboratory    (NSSL)    located   at   Norman, 
Oklahoma,    has    been  gathering   vast   amounts    of   data   during    the    thunder- 
storm season   for    the    past    several   years.      Many   studies    have    been 
conducted  with   the   aid   of   this   data    (for  example   see  Barnes,    1968, 
Newton   and   Fankhauser,    1964   and   Hammond,    1967),    however,    after   an 
extensive    literature    search,    it   appears    that    no   research   has    been 
done,    or   at    least   reported,    involving   a    numerical   modelling   approach 
whereby   the   distributions    of   various    parameters   were   analyzed   and    their 
values    numerically    forecast    for   a    future    time.      Such   an  approach,    if 
successful,    could   give   distinct   advantages    in   predicting,    for   example, 
where   hail   might   be    expected,    or  where    the   most    favorable    conditions 
for    tornadoes   might   be    located.      In  addition,   during   severe    convective 
activity  radiosonde   data    often    terminates    near  mid-tropospheric    levels. 
Consequently,    a   model  which   could    reasonably   predict    the    sounding 
above    termination    level,   when  given  a    complete    sounding    from   the 
previous    observation   time,    would   aid   greatly    in    later   research   efforts. 

Figure    1   shows    the    region   used    for    this    study  with   a   height-field 
analysis    and    the    reported  winds    at    the    nine    radiosonde   stations    in 
the   NSSL  upper-air   network    (see  Appendix  A   for   a  description  of   the 
NSSL  upper-air   network).      Assuming    the   data    to   be    free    of   errors, 
the    non-geostrophic    nature    of   the    flow   on   this    scale    is    immediately 
evident,    and    the    figure    emphasizes    the    need    for    intensive    research 
into   the   dynamic   mechanisms    involved    under   such   circumstances.       In 
an  effort    to   evaluate    the    possibilities    of   numerical  weather    pre- 
diction   for   such   small-scale   data,    this    paper   uses    a    numerical 
prediction  model   based    on   the    complete    isobaric   vorticity   equation. 

13 


Figure    1:      Height    analysis    with   wind   directions    and    speeds  (m   sec 

at    nine    radiosonde    stations    of   NSSL  network    for    the    500-mb  level     it 

1700Z    on    30  May    1967    which    illustrates    the    non-geos trophic  nature    o 
the    data. 
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The  resulting  forecasts  are  then  used  as  input  to  a  cumulus  convection 
model  developed  by  Weinstein  and  Davis   (1967)  to  provide  an  aid  to 
severe  storm  forecasting. 

The  parameterized  cumulus  convection  model  by  Weinstein  and  Davis 
(1967)  calculates  cloud  top,  profiles  of  vertical  velocity,  tempera- 
ture excess  (that  is,  the  temperature  of  a  parcel  of  air  in  the  cloud 
minus  the  temperature  of  a  parcel  of  air  in  the  environment),  hydro- 
meteor  water  content  and  updraft  radius  of  the  cloud  as  well  as  rain- 
fall amount  and  duration.   The  principles  of  the  above  calculations 
are  based  on  the  entrainment  principles  of  Stommel  (1947)  and  the 
results  of  Kessler  (1967).   With  a  standard  radiosonde  sounding  at 
mandatory  and  significant  levels  beginning  at  cloud  base,  the  W.D. 
model  transforms  the  sounding  to  one  which  is  spaced  in  equal  incre- 
ments of  height.   After  the  user  specifies  certain  boundary  conditions 
such  as  ice  nucleation  temperature,  radius  of  the  base  of  the  cloud 
and  conversion  and  collection  rates  (Kessler,  1967),  the  model 
specifies  cloud  extent.   One  use  of  this  model  is  to  compare  the 
effect  of  seeding  a  given  cloud  or  allowing  it  to  develop  naturally. 
This  is  accomplished  by  varying  the  ice  nucleation  temperature. 

The  numerical  model  presented  in  this  paper  predicts  heights, 
temperatures  and  relative  humidities  at  individual  grid  points  for 
seven  isobaric  levels  (950  mb ,  850  mb ,  700  mb ,  500  mb ,  300  mb ,  200 
mb ,  and  150  mb)  so  that  predicted  information  can  be  used  as  input  to 
a  modified  subroutine  form  of  the  W.D.  model  to  yield  the  vertical 


The  model  developed  by  Weinstein  and  Davis  has  been  referred  to  as 
the  Pennsylvania  State  model.   In  this  paper,  it  will  be  abbreviated 
as  the  W.D.  model. 
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profiles  previously  discussed,  under  the  assumption  that  cloud  base 
is  at  850  mb.   However,  the  present  version,  as  described  in  this 
paper,  does  not  provide  for  vertical  coupling  of  the  seven  isobaric 
levels  after  the  first  time  step. 

Predictions  were  made  at  1.5-hour  intervals  utilizing  the  IBM 
OS/360  (MVT  versions  15,  16)  computer.   The  1.5-hour  prediction 
interval  was  selected  to  correspond  to  the  times  when  verification 
data  were  available  from  the  NSSL  radiosonde  network. 

The  grid  used  for  the  computations  is  a  24  x  24  X-Y  grid  oriented 
000-180,  090-270,  having  a  mean  latitude  of  35N.   The  grid  distance 
is  5  nautical  miles  with  the  grid  arranged  so  that  the  nine  upper- 
air  stations  are  either  located  exactly  on  a  grid  point  or  are  so 
close  to  one  that  they  may  be  approximated  to  be  at  the  nearest  grid 
point.   The  farthest  distance  away  from  a  grid  point  for  a  particular 
station  was  2  nm.   They  were  assumed  to  be  located  at  the  nearest 
grid  point  to  eliminate  the  need  for  an  interpolation  scheme  in  the 
objective  analysis  used  for  initialization  of  the  data  as  described 
in  Section  4. 
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2.      MODIFICATIONS    TO   THE  W.D.    MODEL 

While   working  with   and    studying    the   W.D.    model,    it   became   apparent 
that    several  modifications    could   be   made    to    it    that   would    provide   more 
accuracy.      Although   most    of    the   modifications   were    of   a   minor    nature, 
three    of    the   modifications   are    considered    to   be    significant    improve- 
ments   over    the    original  model.      These    three   modifications   will   be 
presented    in  detail  while    the    remainder   of   the    changes   will   be    sum- 
marized  at    the   end    of    this    section   since    they  are    considered    to   be 
self-explanatory. 

Kessler    (1967)    proposed    the    equation 

dM  dm  ,         N  /tin 

dt"         dT  =  ki(m"a)  (2"la) 

to  model  the  conversion  of  cloud  water  to  hydrometeor  water.   Here, 
k1  is  an  empirical  constant  which  is  equal  to  .001  sec   when  m>a 
and  is  zero  for  m<a .   The  quantity  "a"  is  defined  as  a  threshold 
value  and  refers  to  the  cloud  water  content.   If  the  cloud  water 
content  is  greater  than  this  threshold  value,  then  conversion  of 
cloud  water  to  hydrometeor  water  will  occur.   If  cloud  water  amounts 
less  than  the  threshold  value  are  present,  then  conversion  of  cloud 
water  to  hydrometeor  water  will  not  take  place.   This  equation  is 
incorporated  into  the  W.D.  model  (Weinstein  a,nd  Davis,  1967)  in  the 
form 

dQh     dQc 

dr  =  -dr  =  ki<Va'>  <2-lb> 

where  k1  is  defined  in  a  similar  manner  as  it  was  in  (2-la).   Kessler 

(1967)  indicated  that  measurements  have  been  performed  which  indicate 

-3 
that  cloud  water  amounts  greater  than  1  gm  m   are  usually  associated 
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with  the  production  of  precipitation.   Furthermore,  he  arbitrarily 

-3 

suggested  the  value  of  0.5  gm  m    for  the  threshold  value.   Weinstein 

and  Davis  (1967)  chose  a'  =  0.5  gm  Kg  "  for  their  threshold  value. 

It  is  proposed  to  consider  this  threshold  value  to  be  a  function 
of  density.   This  approach  thereby  associates  a  different  threshold 
value  with  each  level  by  making  it  a  function  of  height.   Using  this 
approach,  the  conversion  of  cloud  water  content  to  precipitation  in 
the  W.D.  model  has  been  modified  to 

dQh 

where  p  =  air  density.   By  defining  k...  as  it  was  in  (2-la),  conversion 

of  cloud  water  to  precipitation  is  modelled  to  occur  when  (Q   -  — )  >  0. 

c    p 

Computations  have  been  made  to  determine  the  variation  of  the 
quantity  a/p.   These  values  show  a  maximum  of  approximately  .00137 
gm  gm  "  in  the  upper  levels  and  a  minimum  value  of  .00076  gm  gm 
near  cloud  base. 

The  second  significant  modification  is  the  inclusion  of  a  linear 
variation  of  the  latent  heats  of  vaporization  (L  )  and  fusion  (Lf) 
with  temperature.   Since 

L  =  L  +  Lr  (2-3) 

s    v    f 

and    L     may  be    considered   constant   while    L..    is    a    function   of   temperature, 
s  f 

then   L      is   also  a    function   of   temperature.      The   values    of   L..    (List, 
v  f 

1951)   are   48.6   cal  gm"      at    -50C   and   79.7   cal  gm"      at    0C.      The    latent 
heat   of  vaporization  varies    from  629.3   cal  gm   '    at    -50C   to   597.3 
cal  gm    '    at    0C.      Since    L     varies    by   only    1   cal   gm        from   0C    to    -40C , 
its   value   was    taken  as    constant    and   equal    to   678   cal  gm      .      The 
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variation  of  only  one  of  the  latent  heats  on  the  right  side  of  equation 
(2-3)  needs  to  be  computed  since  the  value  of  the  second  one  can  be 
obtained  by  subtraction. 

In  order  to  approximate  the  variation  of  the  latent  heat  of  vapor- 
ization as  a  linear  function  of  temperature,  it  was  necessary  to  deter- 
mine an  optimum  value  for  the  difference  between  the  specific  heats  of 
water  and  ice.   Integrating  the  formula 

dL 

-PT  =  C    -  C  (2-4) 

dT     pv    w 

and  applying  the  mean  value  theorem  for  integrals  yields 


L  =  L   +  (C    -  C  )  (T-T  ).  (2-5) 

v    v      pv    wy  v   o 
o     r 

Utilizing  the  values  of  L   (List,  1951)  for  various  temperatures  (T) 


and  with  T  =  -40C,  the  value  of  (C    -  C  )  that  satisfies  equation 
o  pv    w 

(2-5)  was  computed  for  each  temperature  at  10  degree  intervals  from 


-40C  to  +40C.   A  plot  of  (C    -  C  )  versus  temperature  is  shown  in 

pv    w  r 

Figure  2.   Based  upon  the  results  it  was  decided  to  select  the  value 


of  -.62340  as  an  "optimum"  value  for  (C    -  C  )  and  consider  it  to  be 

r  pv  w 

a    constant.      This    value    is    considered    to   be   a    representative   value 
between   -40C   and   0C ,    and   corresponds   exactly   to   the   value   at    -IOC. 

The    final  significant   modification   involves    the   values   used    for 
the    constants   when  computing   saturation  vapor   pressure    (Weinstein 
and  Davis,    1967).      Integrating    the   Clausius-Clapeyron  equation  one 
obtains    the    following    for    the   saturation  vapor   pressure    over  water: 


e    [L       -T(C        -   C   )]        e  (C        -   C    )    In  T 

,                                vo           o        pv          w        ,  pv          w 

In  e      =   -  — — L +  — ^ 

s                                     RT  R 


e    [L       -  T      (C        -  C    )]        e    (C        -   C    )  In  T 

vo  o  pv  w  pv  w  o         ,  ._    ,. 

__ — t _  t: — ^ (-   In   e      .  (2-6) 

RT  R  so  ^    u; 

o 
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Simplifying   yields 


where 


lne      =-|~-B'    lnT+C 
s  T 


e    [L     -T   (C        -  C   )] 

,    _  vo      o      pv w 


A'    = 


=  6958.9262 


B'    =   - 


e    (C        -  C    ) 
pv  w 


=  5.65567 


and 


,-L 


c' =1 


vo 


when 
and 


-r^  -    (C        -  C    )    (1+ln  T   )      +   In  e        =  59.01383 
-T  pv  w  o   J  so 

o  r 

lo.  -1 


T     =  233. 16K,    (C        -   C    )   =   -0.62340  cal  gm       K 
o  pv  w 

L        =   621.7    cal  gm 
vo 


Similarly,    when   freezing    occurs    one    obtains 

-A" 


in  %  -  -r +  C" 


(2-7) 


where 


A"   = 


eL 

£ 

R 


=  6151.0205 


and  C"  =  ~-  +  In  e   =  24.3277  when  T  =  233. 16K  and  L  =  678  cal  gm"1. 
T        so  o  s 

o 

The  above  equations  yield  the  saturation  vapor  pressure  in  millibars. 
Obviously  the  number  of  significant  digits  shown  above  exceeds  the 
accuracy  of  some  of  the  assumptions  made  to  arrive  at  the  values  shown. 
Rounding  to  five  significant  figures  could  be  done  with  no  loss  of 
accuracy,  however  they  have  been  retained  here  as  they  were  computed 
on  the  IBM  OS/360  computer. 

There  were  several  minor  changes  made  to  the  W.D.  model.   These 

included  changing  the  value  being  used  for  the  acceleration  of  gravity 

-2 
from  9.87  m  sec    to  a  more  exact  form  obtained  by  computing  the  value 

from  the  formula  (List,  1951)  expressing  g  as  a  function  of  latitude, 

altering  the  values  used  for  the  constants  for  conversion  from  degrees 
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Celsius  to  degrees  Kelvin,  conversion  from  meters  to  feet,  and  the 
value  of  the  gas  constant  for  dry  air.   In  addition,  the  section  of 
the  program  dealing  with  corrections  due  to  the  wind  shear  were  com- 
pletely eliminated  which  consequently  allowed  a  subroutine  (subroutine 
CHEQ)  to  be  eliminated.   Other  minor  changes  included  removing  factors 
of  10  so  that  all  pressure  values  are  carried  in  units  of  millibars, 
introducing  symbols  wherever  possible  to  simplify  future  modification 
and  testing,  changing  the  estimate  of  the  updraft  velocity  at  cloud 
base  from  2.0  to  0.5  m  sec   ,  and  modifying  certain  format  statements 
to  improve  the  headings  on  the  output.   Finally,  many  soundings  were 
run  through  the  final  modified  version  of  the  W.D.  model  to  ensure  that 
all  changes  had  been  properly  programmed  and  to  ensure  that  the  results 
from  the  model  were  reasonable. 
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3.   MATHEMATICAL  DEVELOPMENT  OF  A  FORECAST  MODEL 

This  section  describes  the  model  used  for  this  research.   Although 
it  is  not  consistent  in  many  respects,  it  did  provide  a  means  by  which 
forecasts  of  height  at  1.5-hour  intervals  could  be  obtained  for  input 
to  the  W.D.  model  as  well  as  provide  a  method  by  which  the  values  of 
the  various  terms  in  the  prediction  equation  could  be  ascertained.   A 
description  at  the  end  of  this  section  describes  an  attempt  to  set  up 
the  prediction  equation  in  a  consistent  manner  after  the  original 
model  was  running  properly,  however,  the  more  consistent  method  failed 
to  provide  reasonable  forecasts. 

A.   THE  VORTICITY  EQUATION 

Following  the  procedures  given  by  Haltiner  and  Martin  (1957),  the 
isobaric  vorticity  equation  may  be  written  as 


L(c  +  f)  =  .(C  +  f)v.w  +  (!"!).     <3.U) 


Since  the  entire  grid  used  for  this  research  covers  only  2  degrees 
of  latitude,  the  Coriolis  parameter  can  be  considered  constant, 
taking  on  the  value  it  would  have  at  the  mean  latitude  of  the  grid. 
Using  this  approximation  and  expanding  the  left  side  of  (3-la) yields 

|  .  .  ^|  +  w.vc  +  (c+yv.  w  +  |  |  .  |  |).      (3-  lb) 

In  order  to  use  the  W.D.  model  at  1.5-hour  intervals,  height  of  the 
base  of  the  clouds  is  needed  from  which  computations  can  be  started. 
To  obtain  these  heights,  (3-lb)will  be  used  as  a  height  forecasting 
equation  analogous  to  methods  applied  to  large  scale  motions  in  the 
atmosphere  by  applying  the  balance  equation.  An  alternative  method 
without  using  the  balance  equation  for  initialization  was  tried  and 
is  discussed  at  the  end  of  this  section. 
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Since  the  application  of  (3-lb)  to  height  forecasting  will  be  done 
to  obtain  information  regarding  the  orders  of  magnitude  of  the  various 
terms  and  to  obtain  heights  which  can  be  used  as  cloud  base  information 
for  input  to  the  W.D.  model,  no  attempt  was  made  to  set  up  a  consistent 
set  of  equations.   To  obtain  an  equation  which  could  be  used  for  height 
forecasting  and  evaluation  of  the  various  terms  with  mesoscale  data, 

the  relative  vorticity  on  the  left  side  of  (3-lb)  and  in  the  diver- 

2 
gence  term  was  replaced  by  V  \|f .   In  addition,  only  the  rotational  wind 

components  were  used  for  the  horizontal  advection  term,  that  is, 

u=-|andv=|i.  (3-2) 

This  was  done  because  the  present  version  has  no  provision  for  pre- 
dicting the  divergent  component  of  the  wind  (  \V  ) .   Future  modifications 

X. 

should    provide    the   means    by  which     W      can  be    used    in   addition   to    fore- 

X 

casts  of  the  rotational  component  (  W . ) .   The  vorticity  in  the  vertical 
advection  term  as  well  as  divergence  and  the  u  and  v  components  of  the 
wind  in  the  twisting  term  were  computed  based  on  the  observed  winds 
at  the  initial  time  and  held  constant  for  1.5  hours.   Keeping  values 
constant  in  the  -^—  terms  greatly  simplifies  programming  considerations 
since  3-minute  time  steps  were  used  and  to  update  these  values  at  each 
time  step  would  involve  working  with  each  level  for  only  one  time  step 
before  proceding  to  the  next  level.   Unfortunately,  holding  Q    in  the 
vertical  advection  term  and  the  u  and  v  components  of  the  wind  in  the 
twisting  term  constant  for  1.5  hours  causes  the  seven  levels  to  lack 
vertical  coupling.   Future  modifications  should  attempt  to  correct  this 
deficiency.   After  making  these  substitutions  and  assumptions,  the 
resulting  prediction  equation  becomes 

^2*)="C»|-|I  +  l!|+^+yv.W+^|-||).  (3-3) 
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Assuming   continuity    in   i|j ,    the    term   on    the    left    side    of    (3-3)   may 

9     A I 

be   written  as   V     ■*?*•,    and    the    equation   then  becomes    a    Poisson-type 
equation  which    can   be    solved   by  relaxation   procedures    as    described   by 
Haltiner   and   Martin    (1957).      The    forecasting    procedure,    being   one    of 
successive    iterations,    will    thus    allow  all   terms    involving    ty    in    (3-3) 
to  vary   at   each    time    step. 

To   obtain  a   stream   field    for    (3-3),   consider   the   balance   equation 
(Charney,    1955): 

V*  +  ?  (v^vf)  -  f  L  C^xV  -  ^i  ^fJ  ■  —  •  (3~4) 

Since  the  Coriolis  parameter  is  to  be  considered  constant  over  the 

grid,  the  second  term  of  (3-4)  drops  out. 

As  discussed  by  Charney  (1955)  and  Bolin  (1955),  equation  (3-4) 

will  maintain  its  elliptic  character  only  if 

V  $    f 

^i  +  |>0,  (3-5) 

that  is,  only  if  the  geostrophic  relative  vorticity  is  greater  than 
-f/2.   As  shown  in  Appendix  B,  heights  would  have  to  be  read  at  least 
to  the  nearest  centimeter  on  this  size  grid  in  order  to  satisfy  the 
requirements  of  (3-5)  and  still  allow  negative  geostrophic  relative 
vorticities.   Consequently,  (3-4)  was  forced  to  be  elliptic  by 
writing  the  non- linear  terms  in  the  form 

cM[J\     du  dv  _  ,  . 

^J  =  ~  ^  ^  '  (3"6a) 

and  2    2 

o  f  d  f     dv  du  ,,  ,,  . 

72  72  =  -^^    '  (3"6b) 

dx  dy 

Making  this  substitution  into  (3-4)  and  taking  into  account  that 

the  Coriolis  parameter  is  to  be  considered  a  constant,  (3-4)  becomes 

V2*  =f"  V2z  -  j-  J    (u,  v),  (3-7) 

m        m 
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where    J   represents    the   Jacobian   operator,    and   u  and   v   are    the    components 

of   the    observed  winds.      Equation    (3-7)    is    then  an   elliptic   equation  as 

it   stands   and   does    not    require   that    the  geostrophic   relative   vorticity 

be   greater    than    -f/2.      The    i|f    field   can  be    initialized   by   the   approxi- 

mation   \|r   =  ~— ,    and    then   equation    (3-7)    can  be    relaxed    to   convergence 

m 
before    using    (3-3)    to   predict    the    |    fields    at   a    later    time.      Similarly, 

once    the    ty    fields    have   been   forecast   ahead,    the    height    fields    can  be 

recovered   by  making    the   approximation   z   =  and    then   utilizing    (3-7) 

in   the    relaxation    procedure    until   convergence    is    attained. 

Convergence   was    defined   as    occurring  when    (3-7)    balanced  within 

0.5  meters.      When   converting    from  height    fields    to   stream   fields 

2  -1 

this    requires    an  epsilon   of  approximately    58612   m      sec       ,   whereas 

epsilon  becomes    0.5   m  when   converting    from  stream   fields    to   height 

fields.      These  values   arise    from  the   approximation  \(r   =    (g/f   )z.      In 

retrospect,    the    above    convergence    criteria  may   be    too    large    for    this 

scale    of  motion.      In   the    future,    more    stringent    criteria    (such   as    0.1 

m)    should   be    invoked. 

As   mentioned   at    the   beginning    of   the    section,    attempts   were   made    to 

correct    the    deficiencies    of    the    prediction  equation.      This    involved 

2 
utilizing   the   exact    relation  V  i|/   =  Q   as    the   method   of   obtaining   the 

initial   ty    field    instead    of   using    the   balance    equation  and    to   substitute 

2 

V   y    for   Q    in   the    horizontal   advection   term.      In   addition,    the    conver- 

2  -1 

gence   criteria   was    changed    to    10000  m      sec       ,    in   units    of    f  >   when  con- 
verting   from  (    to   f. 

The    balance    equation  was    used    to   recover   height    fields    from   the 
forecast    fields    of   \|f  with    convergence    criteria    of   0.1   m.      Two  methods 

of   initializing    the    \|f    field   were    tried;    they  were    \[f   =  ~—  and    i|»   =  ~— , 

m  m 
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where  z  represents  the  standard  atmosphere  value  of  the  height  of  the 
respective  isobaric  surfaces.   When  using  z  for  the  initialization, 
computations  were  carried  out  without  difficulty,  however  the  result- 
ing forecast  height  fields  were  very  flat  as  might  have  been  anticipated, 
However,  the  forecast  fields  were  in  error  by  as  much  as  300  m  at  some 
levels.   When  using  z  for  initialization,  instability  resulted  in  the 
horizontal  advection  term  which  caused  an  exponential  growth  of  this 

term  to  values  too  large  for  the  computer  to  handle  after  only  12-15 

2 
time  steps  (36-45  minutes  ahead  in  time).   The  V  if    part  of  the  hori- 
zontal advection  term  is  believed  to  be  the  cause  of  the  instability 
but  this  was  not  investigated  in  detail.   However,  one  attempt  to 
investigate  the  instability  was  made  by  shortening  the  time  step  to 
one  minute.   This  delayed  the  instability  somewhat  which  would  indicate 
the  problems  were  caused  by  computational  instability. 

B.   VERTICAL  MOTION  COMPUTATIONS 

In  order  to  solve  (3-3),  vertical  motion  values  are  required  for 
the  vertical  advection  of  relative  vorticity  and  for  the  twisting  term. 
Vertical  motions  for  each  of  the  seven  levels  were  computed  by  means 
of  the  continuity  equation  in  the  form 

cYjo 


Tf   =  -  V  W  (3-8) 

where  the  bar  symbol  denotes  a  vertical  average  through  the  layer. 
This  method  of  computation  utilizes  the  measured  wind  fields  and  does 
not  allow  the  to  values  to  be  updated  for  forecast  times  since  the 
present  model  has  no  provisions  for  obtaining  the  divergent  component 
of  the  wind.   The  vertical  boundary  conditions  imposed  were 

cd100(J  =  0.0  (3-9a) 

and  cd150  =  0.0  (3-9b) 
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150 
which  also  imply  that  V*  W     =  -  S  V-  \V. 

1000   p  =  950 

Vertical  motion  computations  were  made  by  working  downward  from 

150  mb  through  300  mb  and  upward  from  1000  mb  through  700  mb.   Vertical 

motion  at  500  mb  was  then  computed  as  the  average  of  the  values  at  300 

mb  and  700  mb .   It  is  recognized  that  using  a  backward  or  forward  finite 

cu) 
difference  form  for  -r—  could  at  times  lead  to  computational  difficulties, 

op  r 

however    the   method    used   here   does    allow   the   vertical  motion  at   a   given 
level   to   be   affected   by    not    only   the   divergence   at    that    level  and    the 
level   above,    but   also   takes    into   consideration   the   divergence    field   at 
the    level  below. 

Computations    for    the   various    levels   were   as    follows    (symbols    such 
as  W150   and   W950  were    chosen   to   conform   to   the    symbols    used    in   the 
computer   program    (subroutine   VERTMO)    shown   in  Appendix  C) : 

(1)  150  mb 
W150=0.0 

(2)  200  mb    (Ap>0) 


W175=W150-APV.  W200_150 


W225=Wl50-ApV-  W3()0_150 
W200=(Wl75+W225)/2 
(3)    300  mb    (Ap>0) 


W250=W200-APV.  VV300-200 
W350=W200-APVTW500_200 

W300=(W250+W350)/2 
(4)  1000  mb 

wioor-f .o 
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(5)    950  mb    (Ap<0) 


W925=W1000-APV.  W850_1000 


W975S=W1000-APV.  NV950_10(}0 
W950=(W925+W97  5)/2 
(6)    850  mb    (Ap<0) 


W825=W950-ApV-  W700_950 


W900=W950-APV.  W850_950 
W850=(2 (W825)+W900)/3 
(7)    700  mb    (Ap<0) 


W675=W850-APV.  W500.850 


W725=W850-ApV-  W500_950 
W700=(W67  5+W725)/2 
(8)    500  mb 

W500=(W300+W700)/2 
The   values    obtained    for   vertical   motion   by   this   method   are   discussed 
in  Section   5. 

C.      TEMPERATURE  AND   RELATIVE   HUMIDITY  FORECAST   SCHEME 

In   order    to  utilize    the   W.D.    model   as    an  aid    in   severe    storm   fore- 
casting,   forecasted   values    of   temperature   and   relative    humidity  are 
necessary   in   sounding    form  as    input    to    the   model. 

Although   a   model  designed    for    operational   use    should   definitely 
be    internally   consistent,    such  as    the   vertical   consistency    of    tempera- 
ture   fields   with    the   geopotential   fields,    this    study    involved   deter- 
mining   some    of    the    problems    that   would   develop  when  applying   numerical 
forecasting    techniques    to  mesoscale   data.      Consequently,    the    scheme 
for   obtaining    temperature    forecasts   which    follows    is    not   dynamically 
consistent  with  what   has    preceded,    but    it   does    provide   a   method    of 
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obtaining  a  temperature-field  forecast  which  could  be  used  as  input 
to  the  W.D.  model.  As  discussed  in  Section  6,  future  studies  should 
improve  on  this  method. 

Four  different  methods  of  forecasting  temperature  were  tested  and 
the  results  were  then  compared  to  verifying  data.   The  method  which 
verified  best  was  selected  as  the  forecast  method  for  temperature, 
which  was  simply 

%±  =   -  1/2  W-VT,  (3-10) 

where  the  values  for  the  wind  components  were  taken  as  the  measured 
winds  at  observation  time  and  were  not  altered  during  the  forecast 
interval. 

Similar  testing  of  four  prediction  methods  was  also  done  for 
relative  humidity.   The  results  led  to 

-±   =  -  l/2\VVRh  (3-11) 

as  the  prediction  equation  for  relative  humidity  where  the  wind 
components  were  again  taken  as  the  measured  winds  at  observation 
time  and  were  held  constant  for  the  forecast  interval. 

Temperature  and  relative  humidity  prediction  calculations  by 
equations  (3-10)  and  (3-11)  were  performed  in  subroutine  PROG,  a 
listing  of  which  may  be  found  in  Appendix  C.   Results  of  the  com- 
putations, when  compared  to  verifying  data,  are  found  in  Section  5. 

D.   BOUNDARY  CONDITIONS 

To  initialize  the  i|f  field,  the  approximation  f  ~  gz/f  was  used 
as  discussed  in  Section  3A.  When  the  relaxation  was  performed,  the 
outside  row  or  column  of  the  grid  could  not  be  modified  due  to  finite 
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differencing.      As    a    result,    the    boundary  values    of   i|r  were   held   at    the 
value    obtained   by   the    initialization.      The    prognostic   equation    (3-3) 

was    utilized    from  x   =  3,    22   and   y   =  3,    22.      After    initializing    the 

(ft 

pr1-  field  to  zero,  it  was  then  relaxed  for  the  first  time  step.   This 

relaxation  was  performed  from  x  =  3,  22  and  y  =  3,  22  as  mentioned 
above  and  the  outer  two  rows  or  columns  were  maintained  at  their 
original  values  of  zero.   By  this  means,  actual  forecast  values  of 

|  were  obtained  only  from  x  =  3,  22  and  y  =  3,  22.   After  the  first 

hi 

time  step,  the  previous  relaxed  values  for  ■=rI-  were  used  as  the  initial 

estimate  for  the  new  time  step. 

After  obtaining  the  1.5-hour  forecast  for  the  t  field,  the  con- 
version back  to  z  fields  by  means  of  (3-7)  was  accomplished  by  making 
the  initial  approximation  z  ~  f   f/g  only  for  the  points  where  the  i|f 
field  had  been  forecast.   This  then  leaves  the  outer  two  rows  or 
columns  without  forecasted  height  values.   Before  relaxing  the  z 
field,  however,  the  values  at  these  boundary  points  were  set  by 
performing  a  linear  extrapolation  from  the  fourth  row  or  column 
outward  as  follows : 


3    . 

2    . 

z2    =  z3    - 

0*4   -   z3) 

1    . 

Zl   =  Z2    ' 

(z3    -    z2) 

Once  this  was  accomplished  at  all  of  the  grid  points  in  the  two 
boundary  rows  or  columns,  the  z  field  was  relaxed  by  means  of  the 
balance  equation  from  x  =  2 ,  23  and  y  =  2,  23. 

Horizontal  boundary  values  were  not  a  problem  with  the  vertical 
motion  fields  since  they  were  computed  from  x  =  2,  23  and  y  =  2,  23, 
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As  mentioned  above,  the  prognostic  equation  was  utilized  from  x  =  3,  22 

and  y  =  3,  22  so  that  there  were  computed  values  of  omega  at  x  =  2, 

x  =  23,  y  =  2  and  y  =  23  to  use  for  the  terms  «r—  and  ^r— .   As  stated 

ox     dy 

in  (3-9),  vertical  boundary  conditions  for  od   were  that  the  vertical 
motion  would  vanish  at  150  mb  and  1000  mb ,  which  then  implies  that 
the  divergence  at  1000  mb  is  equal  to  the  negative  of  the  sum  of  the 
divergences  from  950  mb  to  150  mb . 

As  will  be  discussed  in  Section  4,  the  initial  data  fields  of 
height,  temperature,  relative  humidity,  wind  direction  and  wind 
speed  were  analyzed  on  the  computer  by  an  objective  scheme.   The 
analysis  method  included  finite  differencing  methods  such  that  the 
outer  row  or  column  was  not  analyzed.   These  boundary  values  were 
set  to  the  same  value  as  the  next  interior  row  or  column. 

The  same  linear  extrapolation  technique  that  was  used  for  the 
height  fields  was  also  used  for  the  boundary  conditions  of  the  fore- 
casted values  of  temperature  and  relative  humidity,  but  it  was  only 
necessary  to  extrapolate  from  the  second  row  or  column  out  to  the 
first  row  or  column. 

E.   TIME  STEPS 

In   the    forecasting    of    the    new   i|f    fields    by    (3-3),    a    3-minute    time 
step  was    used,    thus    requiring    30   time    steps    to   obtain  a    1.5-hour    fore- 
cast.     The   method    used  was    the    leap-frog    scheme   described   by   Haltiner 
(upbl.),   with   a    forward    time    step    used    for    the    first    step.      To 
eliminate    the    possibility   of   obtaining    two  separate    solutions    by    the 
leap-frog   method,    the    two   solutions   were   averaged   every    10   time    steps 
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(30  minutes)  and  then  one  pass  through  the  averaged  field  was  made 
utilizing  a  Laplacian  smoother  with  a  .04  coefficient  for  smoothing. 
The  leap-frog  forecast  method  was  then  continued. 

The  forecasting  of  temperature  and  relative  humidity  by  the  simple 
advection  schemes  described  in  Section  3C  does  not  involve  any 
iterative  type  computations  since  the  wind  components  were  held 
constant  for  the  forecast  interval.   Consequently,  a  single  90-minute 
forward  time  step  was  used.   However,  future  modifications,  which 
should  include  provisions  for  obtaining  rotational  and  divergent  com- 
ponents of  the  wind,  should  also  update  the  wind  field  with  each  time 
step  so  that  forecast  values  of  the  wind  field  could  be  used  for 
temperature  and  relative  humidity  advection. 
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4.   INITIALIZATION  OF  THE  MODEL 

In  order  to  commence  computations  with  the  previously  discussed 
model,  35  fields  of  data  are  needed  (7  levels  each  of  height,  tempera- 
ture, relative  humidity,  wind  direction  and  wind  speed).   Such  a 
volume  of  data  would  be  difficult  to  obtain  by  manual  means  for 
real-time  operational  use  of  a  model.   Therefore,  a  more  realistic 
approach  is  to  use  an  objective  analysis  scheme  by  the  computer. 

Although  many  complicated  and  very  sophisticated  analysis  schemes 
can  be  and  have  been  developed  (for  example  see  Hughes,  1967),  this 
was  not  the  purpose  of  the  research  reported  in  this  paper.   Con- 
sequently, a  simple  objective  analysis  scheme  was  devised  which  uses 
the  Laplacian  operator  and  an  eight-point  averaging  technique  in 
alternating  steps.   Comparisons  of  computer  analyses  and  hand  analyses 
showed  the  patterns  to  be  very  similar.   The  major  difference  between 
the  two  analysis  methods  was  the  gradient  in  the  vicinity  of  closed 
centers.   The  hand  analyses  tended  to  subjectively  spread  the  gradient 
out  over  a  much  larger  area  whereas  this  was  not  so  with  the  computer 
analyses.   Which  analysis  is  more  nearly  correct  is  merely  a  matter 
of  subjective  judgement.   With  only  nine  data  points  within  a  14,400 
square  mile  area  during  conditions  of  severe  storms,  this  question 
will  probably  not  be  answered  until  the  data  network  becomes  more 
detailed.   Nevertheless,  the  computer  analyses  were  judged  to  yield 
adequate  data  fields  for  testing  of  the  computer  program. 

In  order  to  start  the  analysis  a  first  guess  was  needed  at  each 
grid  point.   This  was  accomplished  by  dividing  the  grid  into  nine 
"regions  of  influence"  (see  Figure  3),  each  region  containing  one 
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Figure  3:   Regions  of  Influence 
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of  the  nine  radiosonde  stations.   The  value  of  a  given  parameter  at 
each  grid  point  within  a  region  was  then  set  equal  to  the  known  value 
at  the  radiosonde  station. 

With  the  initial  guess  completed,  the  analysis  was  started.   The 
first  step  was  to  alter  the  value  at  the  grid  points  (except  for  the 
nine  data  points)  by  setting  the  value  at  that  point  equal  to  the 
average  of  the  surrounding  eight  points  according  to  the  following 
diagram: 

4  3     2 

•  •  • 

5  0     1 


p    Pl  +  P2  +  P3  +  P4  +  P5  +  P6  +  P7  +  P8         (4-1) 
0  ""  8 

After  making  one  pass  through  the  grid  in  this  manner,  the  Laplacian- 

type  scheme  was  used  such  that 

P0=   P0+&    (P1+P3+P5+P7    -   4P0),  (4-2) 

where  the  value  of  6  =0.1  was  arbitrarily  selected.   After  making 
ten  passes  through  the  grid  in  this  manner,  the  eight-point  averaging 
technique  was  again  employed.   The  entire  procedure  was  repeated  in 
such  a  manner  that  the  eight-point  average  was  employed  11  times  while 
the  Laplacian  scheme  was  utilized  100  times.   The  resulting  analysis, 
after  boundary  conditions  were  set,  was  the  analysis  used  in  the  pre- 
diction model.   The  boundary  condition  imposed  was  that  the  value  at 
a  grid  point  on  an  outside  row  or  column  was  equal  to  the  value  at  the 
adjacent  grid  point  on  the  next  row  or  column  toward  the  center  of  the 
grid. 
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Figures  4a  and  5a  show  hand  analyses  for  two  different  levels  and 
Figures  4b  and  5b  show  the  corresponding  computer  objective  analyses. 
A  copy  of  the  program  for  the  objective  analysis  scheme  is  included  in 
Appendix  C. 
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Figure  4a:   Hand  analysis  for  950  mb  at  1700Z  on  30  May  1967 
Contours  every  5  meters. 


38 


Figure  4b:   Computer  objective  analysis  for  95(  mb  at  1700Z  on 
3D  May  1967.   Contours  everv  5  meters. 
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5749 


Figure  5a:   Hand  analysis  for  500  mb  at  1700Z  on  30  May  1967 
Contours  every  10  meters. 
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5749 


Figure  5b:   Computer  objective  analysis  for  500  mb  at  1700Z  on 
30  May  1967.   Contours  every  10  meters. 
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5.   RESULTS  OF  THE  COMPUTATIONS 

Two  sets  of  data  have  been  utilized  with  the  previously  described 
model.   The  times  of  the  two  sets  of  data  are  1700Z,  30  May  1967  and 
2130Z,  30  May  1967. 

In  synoptic-scale  models,  the  terms  in  equation  (3- la)  involving 
vertical  motion  and  divergence  have  often  been  neglected  since  they 
may  be  considered  relatively  unimportant  by  scale  analysis  (Haltiner, 
upbl.)  in  the  large  scale  motions  of  the  atmosphere.   However,  when 
dealing  with  cumulus-scale  data,  this  may  not  be  true.   In  order  to 
evaluate  the  relative  importance  of  the  various  terms  in  the  model 
throughout  the  entire  grid,  computations  of  the  orders  of  magnitude 
of  the  various  individual  terms  of  equation  (3-3)  have  been  performed 
for  both  the  2130Z  and  1700Z  sets  of  observed  data  for  the  seven  levels 
in  the  vertical  described  in  Section  1.   Tables  2  and  3  show  the  mini- 
mum and  maximum  values  of  each  term  based  on  the  observed  data,  not 
on  the  forecast  values.   Although  these  tables  only  show  the  minimum 
and  maximum  values,  they  do  point  out  the  fact  that  the  divergence 
term  is  very  important.   There  were  many  instances  when  all  five 
terms  were  of  the  same  order  of  magnitude,  that  is,  when  each  term 
was  just  as  important  as  each  of  the  other  terms.   Table  4  shows  a 
sample  of  the  values  of  each  of  the  five  terms  as  extracted  from  the 
computer  printouts. 

Out  of  4800  values  tested,  ^—  -^—  and  -^r—  -r—  were  of  the  same  sign 

ox  dp      dy  dp  & 

2584  times,  or  about  54%  of  the  time.   Therefore,  in  the  cases  tested, 
they  were  additive  more  often  than  not  which  is  further  justification 
that  they  should  be  retained  for  small-scale  data. 
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Vertical  motion  computations  yielded  considerably  different  values 
for  the  1700Z  data  and  the  2130Z  data.   Table  5  shows  the  maximum  up- 
ward and  downward  vertical  motions  for  each  isobaric  level  at  each  of 
the  times.   Although  the  values  in  Table  5  may,  at  first,  seem  somewhat 
large,  a  transformation  by  the  formula 

Q)  *  -p  gw  (5-1) 

where  p  is  replaced  by  p/RT  from  the  equation  of  state  to  yield 

^  ~  "  RT  gW'  (5~2) 

shows  that  a  value  of  cd  =  200  mb  hr    upward  corresponds  to  a  value 

of  w  ~  .86  m  sec  "  upward  (at  500  mb)  which  is  certainly  a  reasonable 

value  for  cumulus  convection.   Although  values  this  large  were  only 

present  at  a  few  points,  further  research  is  needed  to  determine  what 

reasonable  values  are  for  this  scale  of  motion.   In  general,  the 

vertical  motion  fields  showed  fairly  strong  upward  vertical  motions 

throughout  the  grid  at  1700Z  and  correlated  fairly  well  with  the 

squall-line  that  was  present  at  that  time.   There  were  a  few  areas 

in  which  the  correlation  was  not  too  good  but  this  may  possibly  be 

due  to  the  fact  that  the  nine  radiosondes  were  not  all  launched 

simultaneously  and  a  difference  of  15  minutes  between  launch  times 

could  possibly  result  in  the  bad  correlation  regions.   Although  the 

radar  pictures  indicated  that  dissipation  of  the  squall-line  had 

begun  at  2130Z  by  the  fact  that  isolated  echoes  were  present  rather 

than  a  large  mass  of  clouds,  the  vertical  motion  fields  did  not  show 

extensive  regions  of  downward  motions  except  at  300  mb .   The  2 130Z 

values  were,  however,  considerably  less  in  intensity  than  they  were 

at  1700Z  which  might  have  been  anticipated  from  the  results  shown  in 

Table  5. 
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The  significance  of  the  much  larger  values  of  vertical  motion  at 
1700Z  as  compared  to  2130Z  was  not  investigated  thoroughly,  however 
larger  values  of  divergence  were  evident  (as  would  be  expected)  at 
1700Z.   This  was  particularly  true  at  200  mb ,  and  this  would  tend  to 
be  propagated  downward  to  300  mb  and  500  mb  by  the  computation  procedure 

The  simple  method  of  forecasting  temperature  resulted  in  values 
which  verified  slightly  better  than  a  persistence  forecast.   The 
average  absolute  error  (not  considering  whether  the  error  was  an 
overestimate  or  an  underestimate)  for  97  persistence  temperature 
forecasts  was  1.8C  while  for  105  50$  temperature  advection  forecasts 
the  average  absolute  error  was  1.4C.   The  reason  for  the  different 
number  of  verifications  is  due  to  the  fact  that  some  of  the  1700Z  and 
2130Z  radiosonde  data  terminated  before   reaching  150  mb .   Verifications 
were  performed  only  at  the  nine  radiosonde  stations,  and  only  at  each 
of  the  seven  discrete  isobaric  levels. 

The  simple  method  of  forecasting  relative  humidity  resulted  in 
values  which  verified  essentially  the  same  as  a  persistence  forecast. 
The  average  absolute  error  for  97  persistence  relative  humidity  fore- 
casts was  18.0$  while  for  105  50$  relative  humidity  advection  forecasts 
the  average  absolute  error  was  18.2$. 

The  results  of  the  height  forecasts  did  not  verify  as  well  as  was 
expected.   This  is  undoubtedly  due  in  part  to  the  fact  that  additional 
refinements  are  necessary  in  the  prediction  method  to  obtain  a  con- 
sistent set  of  equations.   Table  6  compares  the  average  absolute 
errors  of  persistence  and  model  1.5-hour  forecasts  for  both  the  1700Z 
and  the  2130Z  sets  of  data.   As  can  be  seen  in  the  table,  the  model 
did  perform  slightly  better  than  persistence  50$  of  the  time. 


48 


Pi 
o 
Pi 

w 

w 
H 

3 

o 

en 

< 
W 

o 


5 


o 

UO 

1-1 

34.2 

o 

1—1 

CO 

I— I 

CO 

m 

m 

00 

m 

o 
o 

CN 

i— i 
m 

00 

I— 1 
co 

m 

o 

CM 

co 

o 
o 

CO 

CO 
CM 

1—1 

CO 
CM 

r-l 

1—1 

CO 

m 

CO 

o 
o 
m 

O 
CM 

m 

00 

CM 

m 

cr> 

i—i 

o 
o 

CM 

CO 

i— 1 

00 
00 

o 
1—1 

CM 
1^. 

o 

UO 

CO 

CO 

1—1 

o> 

00 

r-l 
C> 

o 
m 
on 

m 

-3- 

1—1 

co 
1—1 

ON 

CM 

1—1 

TYPE  OF 
FORECAST 
(VERIFYING  TIME) 

o 
c 

0) 

4-1    /"N 

CO    tM 

•r-l     O 
CO    CO 
S-i    CO 
CD    i— t 

en  ^-" 

N 

r-l     O 

0)   OO 
T3   CO 
O    r-i 

CD 

g 

4-1    ^-N 
CO    N 

•r-l     O 

CO    O 

r-J      CO 

CD    CM 

Cw    s_^ 

r-l     O 

CD    O 
X>   CO 
O    CM 

0J 

-v 
o 

6 

-0 

d 

CO 

01 

o 

c 

01 
4J 
CO 

•r-l 

CO 

U 

a. 

U 

o 

M-l 


CO 


u 

01 

X! 

> 

M 

0) 

•H 

r— 1 

CD 

-C 

u 

•H 

M-l 

V-i 

o 

CO 

X 

CO 

0 

5-i 

CO 

O 

•r-l 

r-l 

u 

C 

CD 

CD 

> 

CD 

CD 

4-J 

CO 

3 

r-l 

1—1 

O 

1—1 

CO 

CO 

XI 

en 

r-l 

o 

CD 

4-1 

M 

CO 

M 

M 

O 

01 

CO 

E> 

r-l 

< 

^J 

TD 

C 

x> 

CO 

OJ 

N 

r— 1 

O 

JO 

o 

co 

!">» 

H 

i— 1 

49 


As  previously  mentioned,  the  W.D.  model  yields  profiles  of  vertical 
velocity,  hydrometeor  water  content  and  temperature  excess.   In  addition, 
total  rain  expected  from  the  cloud,  duration  of  the  precipitation,  up- 
draft  area,  and  a  profile  of  updraft  radius  are  included  in  the  results 
of  the  computations.   Figure  6  shows  a  sample  of  the  graphical  output 
of  the  W.D.  model.   These  results  are  printed  out  in  graphical  form 
by  the  two  subroutines  (GRPHCL  and  GRAPH)  which  were  included  with  the 
W.D.  model  (Weinstein  and  Davis,  1967).   The  horizontal  axis  (shown  in 
Figure  6  at  1/2  scale)  shows  the  scale  to  be  used  for  the  various  para- 
meters and  the  vertical  axis  on  the  left  shows  the  corresponding 
heights  in  meters,  with  the  bottom  figure  (in  this  case  1642)  being 
the  first  increment  (200  meters)  above  cloud  base.   In  the  actual 
computer  output  of  Figure  6,  the  lines  are  represented  by  letters, 
such  that  the  vertical  velocity  curve  is  a  series  of  the  letter  "W" , 
hydrometeor  water  content  of  the  letter  "Q" ,  temperature  excess  of 
the  letter  "T",  updraft  radius  of  the  letter  "R" ,  and  the  zero  axis 
for  the  temperature  excess  curve  is  represented  by  the  symbol  "$". 

The  W.D.  model  has  a  provision  for  setting  a  constant  value  for 
the  updraft  radius  or  for  allowing  it  to  vary  within  a  cloud.   How- 
ever, for  this  research  the  updraft  radius  was  always  considered  to 
be  constant  and  equal  to  1  km.   The  profile  shown  in  Figure  6  was 
computed  based  upon  the  1830Z, 30  May  1967  forecasted  values  of 
temperature,  relative  humidity  and  height  of  the  cloud  base,  and 
by  assuming  that  cloud  base  was  at  850  mb  for  station  CHK. 

Many  soundings  were  used  as  input  to  the  W.D.  model,  and  in  many 
cases  the  values  of  the  temperature  excess  determined  by  the  model 
seemed  to  be  too  large.   Values  as  high  as  4  or  5  degrees  Celsius 
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occurred  quite  often,  with  an  occasional  value  of  7  or  8  degrees 
Celsius.   This  intuitively  seems  to  be  too  high  but  additional 
observations  are  needed  to  determine  if  such  values  actually  exist 
within  cumulus  clouds.   In  addition,  apparent  discontinuities  in 
temperature  excess,  as  evidenced  near  8000  meters  on  Figure  6,  are 
quite  common  in  the  graphical  output  of  the  model.   These  large  tem- 
perature excess  values  and  temperature  excess  discontinuities  occurred 
not  only  with  forecast  soundings  but  also  with  actual  observed  soundings 
from  the  NSSL  data. 

Figure  6  is  based  upon  an  ice  nucleation  temperature  of  -25C  which 
is  assumed  to  be  representative  of  natural  conditions  in  the  atmos- 
phere.  The  W.D.  model  allows  for  changing  the  ice  nucleation 
temperature  to  a  value  of,  say,  -6C  to  represent  the  conditions  that 
might  occur  if  cloud  seeding  were  to  be  undertaken.   However,  this 
is  presently  not  included  in  the  subroutine  version  of  the  W.D. 
model  (subroutine  CLOUDS)  since  the  reason  for  using  it  as  part  of 
a  prediction  model  was  to  aid  in  severe  storm  forecasting. 
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6.   SUGGESTIONS  FOR  FUTURE  DEVELOPMENTS  AND  IMPROVEMENTS 

The  first  modification  to  the  previously  discussed  model,  which 
became  evident  during  testing  of  the  model,  would  be  to  expand  the 
horizontal  grid  dimensions  to  at  least  32  x  32.   Figure  3  illustrates 
that  four  of  the  nine  grid  points  lie  within  the  first  two  rows  of 
the  boundaries.   This  more  or  less  negates  application  of  the  prognostic 
equations  to  these  grid  points  due  to  finite  differencing.   It  is  esti- 
mated that  this  improvement  would  increase  the  storage  requirements  of 
the  program  by  about  50^,  and  since  the  program  is  already  quite  large, 
modifications  would  have  to  be  made  to  allow  for  the  sharing  of  some 
storage  locations  wherever  possible. 

At  present  the  model  lacks  the  proper  provisions  for  predicting 
the  complete  wind  field  and  for  updating  it  after  each  time  step. 
A  consistent  scheme  for  updating  the  wind  field  should  be  included 
such  that  the  predicted  winds  could  be  recovered  every  time  step. 
In  addition,  this  would  allow  predictions  further  ahead  in  time  as 
well  as  allow  for  updating  the  vertical  motion  fields  with  new 
divergence  fields. 

As  discussed  in  Section  3A ,  the  present  version  of  the  forecast 
model  does  not  provide  for  vertical  coupling  of  the  seven  isobaric 
levels  since  the  vorticity  in  the  vertical  advection  term  and  the  u 
and  v  components  of  the  wind  in  the  twisting  term  were  held  constant 
for  the  1.5-hour  forecast  interval.   To  correct  this  deficiency  would 
involve  performing  one  time  step  at  a  given  level  before  proceding  to 
the  next  level.   Also,  additional  storage  would  be  required  to  save  the 
field  which  •*?—  operates  on  for  one  time  step  back  so  that  vertical 
derivatives  could  be  taken  at  each  level  on  fields  which  all  correspond 
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to  the  same  time.   As  far  as  programming  considerations  are  concerned, 
this  essentially  involves  nesting  the  vertical  level  DO  LOOP  inside 
of  the  time  step  DO  LOOP. 

When  the  model  has  been  modified  to  allow  for  vertical  coupling  of 
the  seven  isobaric  levels,  additional  levels  in  the  vertical  could  be 
included  so  that  vertical  layer  depths  wouid  not  exceed  100  mb .   This 
would  require  perhaps  10  or  11  levels  which  would  increase  the  storage 
requirements  of  the  computer  considerably,  however,  this  would  allow 
a  much  finer  detail  for  the  forecasted  soundings.   In  addition,  it 
should  aid  in  the  delineation  of  the  vertical  velocities  at  each 
leve 1 . 

A  consistent  dynamical  method  should  be  included  for  the  prediction 
of  temperature  and  relative  humidity.   The  present  method  is  more 
empirical  than  dynamical  and  better  verifications  should  be  obtained 
if  a  more  realistic  approach  is  devised.   As  an  alternative,  a  better 
constant  (vice  50%)  might  be  tested,  particularly  in  the  case  of 
relative  humidity  where  the  errors  were  negative  much  more  often 
than  they  were  positive. 

In  subroutine  SETUP  (see  Appendix  C),  cloud  base  is  assumed  to  be 
at  850  mb  for  entry  into  subroutine  CLOUDS  (see  Appendix  C).   A  routine 
could  be  developed  to  be  included  in  SETUP  such  that  the  lifting  con- 
densation level,  mixing  condensation  level  or  convective  condensation 
level  would  be  computed  as  the  cloud  base  before  entering  the  W.D. 
model  (subroutine  CLOUDS).   Comparison  of  results  obtained  for  cloud 
base  at  950  mb  and  850  mb  indicate  that  the  cloud  base  and  the  associ- 
ated temperature  and  relative  humidity  are  very  critical  to  the  results 
obtained  from  the  W.D.  model. 
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Once   a    complete    internally   consistent   model    is    set    up,    calculations 
out    past    1.5   hours    could    be    tried    to   determine   at   what    point    such  a 
model    loses    its    skill.      A   few   runs   were   made    out    to   4.5   hours    and    the 
results   appeared    to    indicate    values    of   height   which  were    too    large,    and 
the    temperature    fields    appeared    to   become    unstable. 

Kessler    (1967)    proposed    the    formula    for    the    fall   speed    of   pre- 
cipitation  to   be    of    the    form 

V  =   5.0939  M'125.  (6-1) 
The    corresponding    formula    in   the  W.D.    model    is 

V  =    15.39   Qh'125.  (6-2) 

-3  -1 

The  units  of  M  are  gm  m   whereas  the  units  of  Q,  are  gm  gm   .   In 

order  to  convert  M  to  the  units  of  Q  ,  it  is  necessary  to  divide  by 

the  density  of  air,  p.   Weinstein  and  Davis  (1967)  have  apparently 

selected  a  representative  value  for  p.   It  is  recommended  that  future 

modifications  of  the  W.D.  model  include  the  density  in  the  computation. 

Since  M  =  pQ  ,  and  p  =  p/RT  from  the  equation  of  state,  it  would  be  a 

simple  change  in  the  programming  to  use  the  following  formula  for  the 

hydrometeor  terminal  velocity: 

/PQhv125 

V  =  5.0939  (jpj)  .  (6-3) 

Since  the  duration  calculation  is  based  on  the  hydrometeor  terminal 
velocity  at  cloud  top,  (6-3)  will  permit  values  of  V  which  not  only 
depend  on  the  hydrometeor  liquid  water  content  but  also  will  depend 
on   the    pressure   and    temperature   at    the    cloud    top. 

Finally,    it    has    been   suggested    to    the   author   by   Dr.    J.    D.    Mahlman 
(personal    communication)    that    the    only   solution   to   overcome    the   many 
difficulties   with    this    scale    of  motion  would    be   a   direct    integration 


55 


of  the  primitive  equations.   It  is  believed  that  no  such  research  has 
been  attempted  for  mesoscale  data  and  this  approach  definitely  warrants 
consideration  for  future  attempts  at  numerical  modelling  of  such  scales 
of  motion. 


56 


7.   SUMMARY 

A  numerical  weather  prediction  model  has  been  described  for  use 
with  cumulus-scale  data.   Together  with  the  W.D.  model,  it  provides 
forecasts  of  height,  temperature,  and  relative  humidity  at  seven 
isobaric  levels  as  well  as  vertical  profiles  of  vertical  velocity, 
temperature  excess  and  hydrcmeteor  liquid  water  content  at  1.5-hour 
intervals.   The  model  was  tested  on  two  sets  of  data  and  yielded  fore- 
casts which  were,  in  general,  statistically  about  the  same  as  a  per- 
sistence forecast.   This  can  probably  be  attributed  in  part  to  the 
fact  that  the  present  version  has  no  provisions  for  coupling  the 
seven  isobaric  levels  nor  for  updating  the  vertical  motion  fields 
after  each  time  step.   Another  factor  is  that  the  temperature  and 
relative  humidity  forecast  schemes  were  more  empirical  than  dynamical, 
and  were  not  consistent  witn  the  vorticity  equation.   This  may  account 
fir  the  lack  of  skill  in  forecasting  these  parameters.   Suggestions 
were  made  to  correct  some  of  these  inconsistencies  in  the  event  that 
further  research  is  conducted  with  the  model.   However,  the  possibility 
of  applying  numerical  forecasting  techniques  to  small-scale  data  does 
appear  to  have  possibi  ities  but  extensive  research  is  needed  in  this 
area  before  useful  forecasts  will  be  obtained. 

The  model  requires  approximately  5.5  minutes  on  the  IBM  OS/360 
(MVT  versions  15 ,-  16)  computer  to  yield  a  1.5  hour  forecast.   An 
additional  14  seconds  per  field  is  required  for  the  objective  analysis 
scheme  described  in  Section  4,  thus  requiring  about  13.5  minutes  for 
the  analysis  of  35  fields  of  data  and  the  prediction  of  the  afore- 
mentioned parameters  for  1.5  hours. 
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APPENDIX  A 

THE  NSSL  RADIOSONDE  NETWORK 

There  were  10  radiosonde  stations  within  the  NSSL  region  considered 
for  this  study  which  are  listed  below.   The  grid-point  numbers  associ- 
ated with  each  station  are  shown  on  Figure  3,  where  the  stations  are 
shown  by  their  call  letters. 


Location 

Call    letters 

Name 

Lat. 

Long. 

SPS 

Sheppard   AFB ,    Texas 

33.97 

98.48 

RIN 

Ring  ling,    Oklahoma 

34.17 

97.58 

LTS 

Altus  ,    Oklahoma 

34.70 

99.33 

FSI 

Fort    Sill,    Oklahoma 

34.65 

98.40 

PVY 

Pauls   Valley,    Oklahoma 

34.70 

97.22 

COR 

Cordell,    Oklahoma 

35.30 

98.97 

CHK 

Chickasha,    Oklahoma 

35.10 

97.97 

OKC 

Oklahoma   City,    Oklahoma 

35.40 

97.60 

TIK 

Tinker  AFB,    Oklahoma 

35.42 

97.38 

WAT 

Wa tonga,    Oklahoma 

35.85 

98.42 

All  of  these  stations,  except  OKC,  were  used  in  this  study.   The 
reason  for  not  including  data  from  OKC  is  that  their  data  corresponds 
to  synoptic  data  times  rather  than  at  1.5-hour  intervals. 
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APPENDIX  B 

ON  THE   ELLIPTICITY  REQUIREMENT   OF   THE   BAIANCE   EQUATION 

The   ellipticity   requirement    of   the   balance    equation    (Charney,    1955 
and   Bolin,    1955)   was    given  by   equation    (3-5)   which    is 

v2f       f 

—  +2>0.  (3-5) 

Bolin  (1955)  stated  that  the  data  at  grid  points  which  did  not  satisfy 
(3-5)  was  modified  so  that  (3-5)  was  satisfied.   This  section  investi- 
gates the  extent  to  which  grid  distance  influences  the  requirement  of 
(3-5). 

Expanding  (3-5)  results  in 

2    o  f 


2         2 

d  f 

m 


z  +  ~  >  0  (B-l) 


where 


2 

V  =  finite  difference  form  of  the  Laplacian  operator, 

m  =  map  factor  ~  1.19, 
g   =  gravity  *=>  9.80, 
d   =  grid  distance, 

z      =  height , 

-4 

and    f     =  mean  Coriolis    parameter  =    .836x10 

m 

Upon  substitution  of  the  above  values  for  each  factor,  (B-l)  reduces  to 

V2z  +  .00087d2  >  0  (B-2) 

where  z  is  in  meters  and  d  is  in  nautical  miles.   Substituting  a 
value  of  d  =  5  nm  yields 

V2z  >  -.02175.  (B-3) 

By  measuring  heights  to  the  nearest  meter,  the  only  way  to  satisfy 

2 
(B-3)  is  to  have  V  z  >  0.   To  satisfy  (B-3)  on  a  grid  with  d  =  5  nm 
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requires  z  to  be  measured  at  least  to  the  nearest  .01  m,  and  this  would 
still  not  guarantee  meeting  the  requirement.   It  is  easily  shown  with 

equation  (B-2)  that  the  smallest  grid  distance  on  which  (3-5)  can  be 

2 

satisfied,    without   eliminating    the    possibility   of  V  z   <   0  and    still 

measure    z    to   the    nearest   meter,    is    35   nm. 

The    result    from   above    is  what    prompted    the    substitution   of   the 
measured   wind    field    into    the    balance    equation   in   order    to    force    it 
to   be    elliptic.      Tests   were   made    to   determine   with  what    frequency   the 
balance   equation  ellipticity   requirement   was    satisfied   with   this    scale 
of   data.      The   results    show   that    out    of   6776    points,    only   62.4%   satis- 
fied   the    requirement.      At    1700Z    the    requirement   was    satisfied    59.4% 
of   the    time   while   at    2130Z    the    figure   was    65.5%. 
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APPENDIX  C 


This  appendix  includes  a  listing  of  the  computer  program  used  for 
the  research  described  by  the  preceding  sections.   In  order  to  aid 
any  possible  research  in  the  future  with  this  program,  each  subroutine 
will  be  listed  with  a  brief  explanation  of  the  calling  arguments. 
1.   Subroutine  VERTMO 

U  =  Zonal  component  of  the  horizontal  wind  in  m  sec 

V  =  Meridional  component  of  the  horizontal  wind  in  m  sec 
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OMEGA  =  cd  =  Vertical  component  of  the  wind  in  (x,  y,  p,  t) 

coordinate  system  in  mb  sec 
NI  =  Number  of  grid  points  along  x-axis . 
NJ  =  Number  of  grid  points  along  y-axis. 
NK  =  Number  of  levels  in  the  vertical. 
DIST  =  Distance  between  grid  points  in  meters . 
DSI  =  The  negative  of  the  divergence  field  for  1000  mb  =  the 

sum  of  the  divergence  from  950  mb  to  150  mb . 
Subroutine  FORFUN 

VORT  =  Vertical  component  of  relative  vorticity. 
F  =  Field  for  forcing  function. 
FMAP  =  Map  factor  (mean). 
FBAR  =  Coriolis  parameter  (mean). 
DIST  =  Distance  between  grid  points  in  meters. 
DT  =  Forward  time  step  in  seconds. 

USPD  =  Zonal  component  of  the  horizontal  wind  in  m  sec 
VSPD  =  Meridional  component  of  the  horizontal  wind  in  m  sec 
K  =  The  level  the  main  part  of  the  program  is  at  when  this 
subroutine  is  called. 
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IZ  =  Counter  to  determine  at  what  time  step  (1-30)  the  main 

program  is . 

2     -1 
PSI  =  The  l|f  field  in  m   sec 

OMEGA  =  oj  =  Vertical  component  of  the  wind  in  (x ,  y,  p,  t) 

coordinate  system  in  mb  sec 

3.   Subroutine  RELAXI 

M  =  Number  from  2,  NI  =  a  number  one  greater  than  what  is 

needed  in  DO  LOOP. 

N  =  Same  as  M  except  2,  NJ. 

F   =  Aty    field    for   relaxation   of   balance    equation. 

2  -1 

A  =   ty    field    in  m      sec 

B  =  Z  field  in  meters. 

L  =  Counter  to  keep  track  of  number  of  points  which  have 

converged . 

EPS  =  Epsilon  =  convergence  criteria. 

ALPHA   =   Over-relaxation  constant. 

FMAP   =  Map    factor    (mean). 

FBAR  -  Coriolis  parameter  (mean). 

DIST  =  Distance  between  grid  points  in  meters. 

IL  =  Lower  point  for  DO  LOOP  in  x  direction. 

JL  =  Lower  point  for  DO  LOOP  in  y  direction. 

LM  =  Iteration  number  from  main  program. 

U  =  Zonal  component  of  the  horizontal  wind  in  m  sec 

V  =  Meridional  component  of  the  horizontal  wind  in  m  sec 

K  =  Level  in  the  vertical  from  main  program. 

G  =  Gravity. 
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4.  Subroutine  RELAX 

Same  as  in  REIAXI  except: 
F  =  Forcing  function. 
R  =  Residual  field. 
ICOUNT  =  L. 
II  =  IL. 
Jl  =  JL. 

5.  Subroutine  SMOOTH 

D  =  Field  to  be  smoothed. 

RF  =  Smoothing  factor  (<1.0). 

IPASS  =  Number  of  smoothing  passes  to  make  through  the  grid 

IL,  JL  =  Lower  point  for  DO  LOOP. 

IH,  JH  =  Upper  point  for  DO  LOOP. 

6.  Subroutine  PROG 
TEMP  =  Temperature. 

RH  =  Relative  humidity. 
DIR  =  Wind  directions. 
SPD  =  Wind  speeds. 


USPD   =  Zonal   component    of   SPD   in   m   sec 
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VSPD   =  Meridional   component    of   SPD    in  m   sec 

DIST   =  Distance    between  grid    points    in   meters 

NK   =  Number    of    levels    in   the   vertical. 

NI    =  Number    of  grid    points    on   x-axis . 

NJ   -  Number    of   grid    points    on   y-axis . 

Nil   =  NI-1. 

NJ1   =  NJ-1. 

OMEGA   =  a). 

Z  =  Height  field . 
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7.  Subroutine  SETUP 
TEMP  =  Temperature. 

RH  =  Relative  humidity. 

Z  =  Height  field. 

IL  =  A  number  from  1  to  NI ,  usually  1. 

IH  =  A  number  from  1  to  NI ,  usually  NI . 

JL  =  A  number  from  1  to  NJ,  usually  1. 

JH  =  A  number  from  1  to  NJ,  usually  NJ. 

8.  Subroutine  CLOUDS. 

This  is  the  modified  version  of  the  W.D.  model  (Weinstein  and 
Davis,  1967). 

PDATA  =  Pressure  levels. 

ZDATA  =  Height  field. 

RHDATA  =  Relative  humidity  field. 

ZZ1  =  Height  of  cloud  base. 

DZ  =  Vertical  increment  (meters)  to  be  used. 

NCALL  =  Number  of  times  this  subroutine  called.   Used  as 

variable  to  bypass  a  READ  statement. 
IPT  =  I  coordinate  for  set  of  data  coming  in. 
JPT  =  J  coordinate  for  set  of  data  coming  in. 

9.  Calling  arguments  of  GRPHCL  and  GRAPH  are  all  parameters  which 
have  been  computed  or  set  in  subroutine  CLOUDS. 

10.  Subroutine  METMAP 

This  subroutine  was  not  written  by  the  author,  but  is  part 
of  the  Naval  Postgraduate  School  Computer  Facility  Library  and  has 
been  included  for  information.   It  is  simply  a  shading  routine  that 
will  print  out  a  field  of  data  with  a  0.5  inch  grid  spacing  and  contour 
the  field. 
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Y  =  Two-dimensional  field  to  be  contoured. 

N  =  Number  of  rows  I  in  the  array  to  be  contoured. 

M  =  Number  of  columns  J  in  the  array  to  be  contoured. 

T  =  Title  for  printout  .   Up  to  96  columns  of  alpha-numeric 

information. 
BND  =  Bandwidth  desired  for  contouring. 
AZ  =  Scaling  constant;  that  is,  each  value  will  be  multiplied 

by  AZ  before  printing  it  out.   Only  the  first  three  numbers 
to  the  right  of  the  decimal  point  are  on  the  output,  so  AZ 
can  be  used  to  control  which  numbers  are  to  be  printed  out. 
BZ  =  Additive  constant.   This  is  useful  when  working  with  D-values. 
Usually,  however,  it  is  0.0.   If  used,  this  constant  will  be 
added  to  each  value  of  Y  before  printing  it  out. 
AMIN  =  Minimum  value  for  subroutine  to  start  contouring. 
IJT  =  0  means  subroutine  will  compute  the  minimum  value  and 
start  contouring  at  that  value.   A  value  must  still  be 
specified  for  AMIN  but  it  will  not  be  used. 
ICON  =  1  if  contouring  desired. 

=  0  if  no  contouring  desired. 
When  working  with  this  subroutine,  it  must  be  remembered  that  the 
grid  point  I  =  1,  J  =  1,  is  the  upper  left  corner  of  the  grid,  and  I 
increases  downward  while  J  increases  to  the  right.   If  an  array  is 
defined  with  the  lower  left  corner  as  the  point  (1,  1),  the  follow- 
ing set  of  FORTRAN  statements  will  get  the  array  in  the  proper  order 
for  the  subroutine: 
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DO  5  1=1,  24 
DO  5  J=l,  24 
5  DUMMY  (J,I)=Z(I,24-J+1) 
The  field  called  DUMMY  is  then  taken  into  METMAP  as  Y. 
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